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I. INTRODUCTION 

Thanks to considerable progress achieved in variational 
approaches [1-3], Hamiltonian investigations of Yang- 
Mills theories in Coulomb gauge [4] have seen a renewed 
interest in the last years. In such setup Gaufi's law on the 
states is cxplicitcly resolved with the help of the gauge 
constraint, without any need to construct the physical 
Hilbert space explicitly [5] . The results obtained [6] agree 
with the Gribov-Zwanziger confinement scenario [7-9] 
and with the existence of a confining potential driven by 
topological excitations [10-12]. Lattice calculations have 
confirmed all qualitative and most quantitative features 
of the continuum analysis [13-16]. 

Recently, an extension of variational techniques to full 
QCD has been proposed in Ref. [17]. In the present paper 
we wish, via direct lattice calculations, to set the bench- 
mark for future checks of such results and put forward 
further predictions. 

A key issue in all lattice studies of pure Yang-Mills the- 
ories in Coulomb gauge has been the renormalizability of 
Green's functions. Indeed it has been shown in Refs. [13- 
16] that for static correlators a non-perturbative renor- 
malization procedure can be defined in the lattice Hamil- 
tonian limit at — >■ 0. Such limit must be either taken ex- 
plicitly [16] by going to anisotropic lattices [18], or it can 
be circumvented whenever a decoupling of the temporal 
and spatial dependence of the correlators is possible [13- 
15]. Lacking a clear picture from perturbation theory 
[19-21], this is the only setup we are aware of in which 
multiplicative renormalizability of Coulomb gauge corre- 
lators seems to be guaranteed. 

We will show in this paper that, as for the gluon, the 
lattice Coulomb gauge static quark propagator S{p) = 
J dpiS{p) is indeed renormalizable. From S{p) the rcnor- 
malization function and the running mass M(|p|) 
can be extracted. The mass function Af(|p|), in par- 
ticular, encodes the relevant information on chiral mass 



and the anomalous dimension; we will find that it ba- 
sically coincides, for small momenta, with its Landau 
gauge counterpart. Furthermore, we also demonstrate 
that the full Coulomb propagator S{p) has a trivial en- 
ergy dependence. Besides making it renormalizable, this 
allows for a definition of a quark effective energy com- 
patible with the confining properties of the theory. 

The paper is organized as follows: Sec. II describes the 
general lattice setup used in our investigation; Sec. Ill 
defines the observables we will study and the require- 
ments for their renormalizability; in Sec. IV we will give 
our results and discuss their consequences, while Sec. V 
contains our conclusions and an outlook for future work. 



II. LATTICE SETUP 

A. Gauge field configurations 

Our calculations have been performed on nine sets of 
gauge field configurations generated by the MILC col- 
laboration [22, 23], made available via the Gauge Con- 
nection [24] and analyzed using mostly the C++ toolkit 
FermiQCD [25]. The configurations were produced with 
the Symanzik-improved Liischer-Weisz gauge action [26] . 
Seven out of nine sets include two light degenerate [u, 
d) and one heavier (s) quark flavor, while two are in the 
quenched approximation; all dynamical calculations used 
the Asqtad improved action [27]. The parameters of all 
sets are summarized in Table I; for the reported lattice 
scales and quark masses we refer to Refs. [22, 23]. The 
lower part of the table also contains three sets, (j)-(k), 
of quenched configurations on smaller lattices used for 
cross checks; they were produced with the Wilson action 
at /3 = 5.7, 6.0 and 6.5. For set (i), due to the I/O han- 
dling of FermiQCD, our main bottleneck, we could only 
analyze five configurations in reasonable time; although 
averaging over time slices did considerably increase the 
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statistics, auto correlations still remained high. This is 
reflected in the relatively large error bars on final quan- 
tities for this lattice; we will further comment on this in 
Sec. IV C. 

Throughout the paper a will denote the lattice spacing, 
X and p refer to the lattice 4-coordinate and momentum, 
respectively, while x, p denote the spatial and X4, the 
temporal position and momentum component. Repeated 
Greek indices indicate summation over all four Euclidean 
components while repeated Latin indices refer to spatial 
components only. 





X T 


[fm] 


am 


m [MeV] 


# config. 


(a) 20^ 


X 64 


0.121 


0.010, 0.050 


15.7, 78.9 


202 


(6) 20=^ 


X 64 


0.121 


0.020, 0.050 


31.5, 78.9 


50 


(c) 20^ 


X 64 


0.120 


0.030, 0.050 


47.3, 78.9 


25 


(d) 20^ 


X 64 


0.119 


0.040, 0.050 


63.1, 78.9 


25 


(e) 20^ 


X 64 


0.121 






66 


(/) 28« 


X 96 


0.086 


0.0062, 0.031 


14.0, 67.8 


50 


(S) 28^ 


X 96 


0.086 


0.0124, 0.031 


27.1, 67.8 


50 


(h) 28^ 


X 96 


0.086 






77 


(i) 


X 144 0.060 0.0036, 0.0108 11.8, 35.3 


5 


(i) 16^ 


X 32 


0.170 






94 


(fe) 16^ 


X 32 


0.093 






108 


(/) 16^ 


X 32 


0.045 






131 



TABLE L The gauge field configurations used in this study. 
Sets (a)-(d) and include two degenerate (m, d) quarks 

and one heavier (s) flavour [22, 23]. 



B. Gauge fixing 

The continuum Coulomb gauge condition, 

V • A = (1) 

can be realized on the lattice by maximizing the gauge 
functional 



Fg[U]=mcJ2^r[Uf{x) + U^{x-i) 



(2) 



with respect to gauge transformations g{x) G SU(3), 
where 



U^^{x) = g{x)U^{x)g{x + il)l 



(3) 



The gauge links maximizing Eq. (2) will satisfy the dis- 
cretized Coulomb gauge condition 



where we define the gauge fields as: 
"[/,(x)-C/,(a;)t 



Ai{x) 



2iago 



(4) 



(5) 



This corresponds to Hermitian generators T" = A°/2 for 
the Lie algebra of SU{3), where A" are the Gell-Mann 
matrices. 

To maximize Fg[U] we employ the standard over- 
relaxation algorithm [28], which in the case at hand can 
be applied to each time-slice independently. The inver- 
sion of the Dirac operator still needs however to be per- 
formed on the whole lattice and is computationally quite 
expensive. This forces us to analyze less configurations 
than those used, say, to calculate the Coulomb opera- 
tor [16], making statistics lower and raising the ratio of 
statistical noise to Gribov noise; we have therefore not 
seen, at least for this first study, any need to adopt the 
improved techniques first developed for Landau gauge in 
Rcfs. [29, 30] and adapted to Coulomb gauge in Rcfs. [13 - 
16] . Due to the presence of fermions the part of such tech- 
niques based on center transformations and thus topolog- 
ically non-trivial Zjv sectors [29, 30] could have not been 
applied anyway. 

A measure of the quality of the gauge fixing for each 
fixed time-slice is the average L2-norm of the gauge fixing 
violation ^ [31] 



— y 



tr {ISP{x,x^)N>(x,x^)'^ 



where the sum runs over all spatial sites x and 1? is the 
mimber of lattice sites in one time-slice. We have chosen 
to stop the over-relaxation algorithm whenever on each 
time-slice \fd^ the default output of the FenniQCD code, 
reached machine single-precision, < 510-^ Notice 
that in the literature (see e.g. Ref. [32]) the value of Q 
itself, rather than of is often reported for the quality 
of the gauge fixing; in our case this corresponds to < 
2.5 10"^'^. A possible improvement we are considering to 
implement in the future is to introduce a local rather than 
a global observable in triggering the stopping criterion, 
see e.g. Ref. [29, 30], resulting in a better gauge fixing. 

As usual in Coulomb gauge maximizing Eq. (2) leaves 
the temporal links lJ^{x^ unfixed, i.e. we still have a gauge 
freedom with respect to space independent gauge trans- 
formations g(x\) e SU(3), which leave Fg\U\ unaffected. 
One possible choice to fix the residual gauge in the con- 
tinuum is to require 



^4 ^ ^xAi{x) =0. 



(7) 



traceless 



Throughout this paper we will use the lattice version of 
Eq. (7) proposed in Ref. [13, 14], which we here call inte- 
grated Polyakov gauge (IPG) and which we have adapted 
to the colour group SU (N) . For sufficiently fine temporal 
lattice spacing the IPG might be considered a good ap- 
proximation, on average, to the Weyl gauge ^4 = used 
for the Hamiltonian approach in the continuum. For fur- 
ther details see Appendix A. 
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III. THE QUARK PROPAGATOR 

A. Quark propagator in Coulomb gauge 

At tree-level the inverse continuum quark propagator 
in Euclidean space reads: 

S^''\p)-^ =ii>+ii>i + m, (8) 

where we have used Feynman's slash notation, ;i = 
YliiliPi S'lid f)i = 74^4 ; m denotes the bare quark mass. 
We have explicitly separated the spatial momenta pi from 
the temporal component p4 (the energy) to make contact 
with the non-manifestly Euclidean invariant interacting 
Coulomb gauge propagator S~^{p). The latter could, in 
principle, have a more complex Dirac structure than in 
covariant gauges since three-dimensional (spatial) covari- 
ance allows for more possible contractions of the momen- 
tTim components with elements of the Euclidean Clifford 
algebra. More precisely, the spatial momentum p can be 
contracted with both the spatial Dirac matrices 7^ and 
the proper vector component 0-44 of the Euclidean tensor 
(^nf = V2[7^,7^]. Thus, S^^{p) can be decomposed into 
four pieces: 

S~Hp) = iMsi\p\,P4) + ii>4At{\p\,P4) 

+ 'J2^'^'^:i'^^d{\p\,P4) + Bm{\p\,P4) (9) 
j 

with scalar functions As{\p\,p4), At{\p\,p4), Ad.{\p\,P4) 
and Bm{\p\,P4), to which we will refer as the spatial, 
temporal, mixed and massive component, respectively. 

For our caknilations of the valence quark propagator 
we have used Asqtad improved [27] staggered fermions 
[33], concentrating on the dynamical point for the un- 
quenched configurations. In Landau gauge this choice 
[32, 34-38] showed not only good agreement with the 
conceptually cleaner but very expensive overlap fermions 
[39 45] but also suffered less from short-distance cutoff 
effects as compared to clover (improved Wilson) [46, 47] 
or chirally improved (CI) fermions [48]. 

Only for the finest lattice (set (i)) we were forced to use 
standard Kogut-Susskind fermions [33] due to memory 
limitations. For this set we have calculated the propa- 
gator with a mass equal to the two light degenerate dy- 
namical quarks (11.8 MeV) plus four partially quenched 
masses up to 142.1 MeV. On the quenched configuration 
{h) we have used a mass of 14.0 MeV, i.e. the same as on 
configuration (/), in order to study the effects of dynam- 
ical quarks. 

Notice that for standard staggered fermions S{p) will 
actually be a function of fc^ = sinp^, while the equivalent 
expression for the Asqtad improved quarks can be found 
in Appendix B. To keep notations simple we will how- 
ever write the structure functions Ag, At, A4 and Bm as 
functions of the the discrete momenta p^ taking values 
in the first Brillouin zone. 



B. Dispersion relation and renormalizability 

From the configurations at our disposal, once any suit- 
able gauge has been fixed, the inversion of the Dirac oper- 
ator directly provides us with the regularized propagator 
5reg(a;p), which depends on the lattice spacing a. As- 
suming multiplicative renormalizability, such regiilarized 
propagator should be related to the renormalized one 
Sq{p) via the quark wavefunction renormalization con- 
stant Z2, which will depend on a and the renormalization 
point C, 

5reg(a;p) = ^2(C;«)^c(p)- (10) 

In Coulomb gauge the static propagator can be then 
extracted from Sq{p) by integrating it over p4. For 
bosonic fields such static propagator agrees, up to a con- 
stant, with the inverse of the boson's dispersion relation 
ujb{\p\)- This connection is essential in showing e.g. that 
the Gribov-Zwanziger confinement scenario for the gluon 
is indeed realized within pure Yang-Mills theories [13-15]. 

In the fermionic case things are a bit more complicated. 
In Landau gauge, up to well understood discretization 
effects, the renormalized propagator Sc^{p) was shown [32, 
34-38, 48] to have indeed the expected form 

where the only dependence on the renormalization scale 

C is through the numerator Z^ip^), while the mass fimc- 
tion M{p'^) in the denominator does not depend on the 
cutoff a or the renormalization scale C,; it is thus a renor- 
malization group invariant. 

What is the physical meaning of such propagator and 
where can the dispersion relation be read from? At least 
for the non interacting case, where M = m, and Z is 
constant, we know that the inverse propagator Eq. (11) 
is proportional, up to a Wick rotation, to the projectors 
on the positive and negative energy states [49]: 

K±{p)oc.±i)+m; (12) 

the free particle dispersion relation can now be obtained 
by integrating the square of such Euclidean propagator 
in the energy p4, 

On the other hand, Eq. (11) in the interacting case can 
be considered the propagator of a single quasi-particle, 
where M depends on the 4-momentum p. In this inter- 
pretation a projector as in Eq. (12) can still be defined by 
simply substituting m M{p'^). To extract the effective 
energy of the quasi-particle, however, one needs to model 
the functional form of M and Z in order to perform the 
integral as in Eq. (13), since a numerical summation over 
P4 would be plagued by cut-off effects [13, 14]. 
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In Coulomb gauge, wc can study both the full energy- 
dependent propagator as in Eq. (9) and the static propa- 
gator S{p) = J dpiS{p), taking the form (see Sec. IV C): 



Sdp) 



zd\p\ 



ii>+M{\p\) ■ 



(14) 



What is the meaning of these quantities? Consider the 
free Dirac equation: 

ig^i} = Pitp = h{p)ip = {■}/)+ 7X1)11) , (15) 

where the one-particle Hamiltonian h{p) coincides, up to 
a gamma matrix and a Wick rotation, with the free in- 
verse static propagator (i.e. M = m and Z constant). 
We can again define the projectors A± on the posi- 
tive/negative free energy states and find upon Wick ro- 
tation: 



~A+{p) - A^ip) ^ hip) S-\p) 
A+(p) + A_(p) = l, 



(16) 



while the dispersion relation is given by the eigenvalues 
of h{p). If we now turn on the interactions and appeal 
again to the quasi-particle description we will have in 
general: 



h{p) 



1 



P{\P\ 



-74(/i+M(|p|)) 



(17) 



and the (Wick rotated) static Coulomb propagator will 
still define the energy projectors. Moreover, if we could 
show for the full Euclidean Coulomb propagator a rela- 
tion of the form: 



sdp) = ■ 



zd\p\: 



ip+ipia{\p\) + M{\p\ 



(18) 



we could, up to a constant, directly read from the integra- 
tion of S'^ in p4 the dispersion relation i.e. the effective 
energy of the quasi-particle: 



^ J 27r a^{\p\)pi+p- 



Z!i\p\) 
a{\p\)^p^ + M^{\p\)' 



'2 + M2(|p|) 

(19) 



To check the above conjectures we first need to extract 
the four dressing functions in Eq. (9) from the lattice cal- 
culation. This can be done by extending the techniques 
developed for Wilson and staggered fermions in Landau 
gauge in Rcfs. [34, 37]. For sake of readability all techni- 
cal details have been deferred to Appendix B, where the 
staggered case is explicitely discussed. The adaptation of 
the method to Wilson and Asqtad fermions is straight- 
forward. 



C. Data cuts and averaging 

Throughout this work we average over the cubic (spa- 
tial) symmetries of the lattice and the parity symmetry 
of the propagator S{p) = S{—p). To minimize discretiza- 
tion artifacts caused by the breaking of rotational invari- 
ance we perform a cylinder cut on the spatial momenta 
[50, 51]. 



IV. RESULTS 

A. Dirac structure 

Fig. 1 shows the structure function for some of the 
data set in Tab. I, plus one of the quenched Wilson simu- 
lation, set (fc) (Fig. Id), used as a cross check. In all cases 
Ad = 0. whic;li extends the one loop result of Ref. [52] to 
the non perturbative regime. We will thus ignore the 
mixed component in the rest of this paper and express 
the lattice propagator as: 

S~^{p) = illkaAs{\p\,p4) + ift4aAt{\p\,p4) + Bm{\p\,P4) , 

(20) 

where we denote with fc^ the dimensionless lattice mo- 
menta (sec Appendix B) while the spacing a in the first 
two terms renders the dressing functions Ag and At di- 
mensionless. 



B. Energy dependence 

In Refs. [13 16] it was shown that the static propaga- 
tors of pure lattice gauge theory in Coulomb gauge are 
subject to cutoff effects in the temporal lattice spacing 
at, which lead to violations of multiplicative renormaliz- 
ability. These effects can be both direct, i.e. caused by an 
explicit energy dependence of the correlator at hand [IS- 
IS], and indirect, i.e. caused by 0{at) corrections to the 
spectrum of the theory [16]. In principle we cannot ex- 
clude the latter effect for the quark propagator. However, 
the indirect energy dependence is usually much smaller 
than the direct effect, and a higher statistical precision as 
that in the present work, combined with simulations on 
anisotropic lattices [18], would be required to resolve it. 
As for the direct energy dependence of the quark dressing 
functions, a quantitative measure is given by: 



dm{z) 

dt{z) 
ds{z) 



Bm{\p\,PA) 



At( 



I Pi ) 

,P4) 



A(iPi,pr") 

As{\p\,Pi) 

As{\p\,prn 



(21) 



which on dimensional grounds should only be functions 
of = . In Fig. 2 we show as an example the functions 

of Eq. (21) for configurations sets (a) and (i), plotted as 
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0.25 
0.2 
0.15 
% 0.1 
< 0.05 


-0.05 
-0.1 
-0.15 



Coulomb gauge 
Coulomb gauge + IPG 




0.5 1 



1.5 2 2.5 
|k| [GeV] 



(a) Configs. (a) 



3.5 4 



0.1 



< 



-0.1 



Coulomb gauge 
Coulomb gauge + IPG 



f f 



0.5 1 1.5 2 2.5 3 3.5 4 
|k| [GeV] 

(b) Configs. (a), p4 averaged 



I I 
I 



''lk| (Ge\^) 




-0.2 



-0.4 ■ 



6 a 

p [QeV] 

(d) Configs. (fc), p4 averaged 



10 12 



(c) Configs. (6), p4 averaged 
FIG. 1. Dressing function Ad. In Subfigs. la, lb the effect of residual gauge fixing (IPG) is shown. 



a function of 1 + to better compare with Fig. 1 of 
Ref. [13]. 

We can conclude that, at least within our numerical 
precision, the functions As, At and Bm are independent 
of the energy p4. Contrary to the gluon propagator [IS- 
IS], explicit violations of renormalizability due to the 
temporal cutoff at can therefore be ruled out. If we do 
keep p4 as an argument it is only to distinguish the p4-un- 
averaged structure functions obtained after the temporal 
gauge has been fixed from the p4-averaged ones which 
describe the static functions obtained for a fixed time 
slice. 



C. Full vs. static propagator and their 
renormalization 

Due to the energy independence of the dressing func- 
tions Eq. (20) we can average them over to minimize 
statistical fluctuations. The full propagator thus reads 

S-\p) = i^taAsilpl) + ihaAti\p\) + B„A\p\) , (22) 

while its static counterpart is, up to a constant propor- 
tional to the time extent of the lattice (see Appendix C): 

S-\p) ^ iltaAs{\p\) + BM) ■ (23) 



Alternatively, the static propagator can be taken at a 
fixed time slice, which is equivalent to averaging over 
all time slices without residual gauge fixing. As shown 
in Fig. 5a, At identically vanishes in this case and we 
again obtain Eq. (23). Such fixed time definition would 
of course make the static quark propagator much eas- 
ier to calculate on lattices with large temporal extension, 
since for xT the inversion of the Dirac operator can be 
restricted to the spatial sub-lattice. Even if the proce- 
dure is repeated T times to improve the statistics one will 
still cut both CPU time and memory requirements and 
avoid possible I/O inefficiencies. We are planning to im- 
plement such improvement in all future analysis; besides 
enabling us to increase the number of configurations, this 
should also allow us to invest more computer time in the 
quality of the gauge fixing, see discussion in Sec. II B. 
Of course, this does not apply to the full propagator of 
Eq. (22), which will still need the inversion of the whole 
Dirac operator to resolve Af. 

Let us now define the renormalized propagators. While 
from Eq. (23) a static propagator as in Eq. (14) imme- 
diately follows, an explicit p^ dependence of the dressing 
functions would have made it very difficult to cast S{p) 
to the form Eq. (18), with the scale dependence con- 
fined to the proportionality factor. Due to the trivial 
P4-dependence of Eq. (22), however, the full (non-static) 
quark propagator does indeed take the lattice equivalent 
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(a) dm from (a) 



(b) dt from (a) 





(c) from (a) (d) ds from (i) 

FIG. 2. Energy dependence of dressing functions as defined in Eq. (21). 



form of Eq. (18): 

sdp) = 



zd\p\ 



(24) 



iafk + iajtiailpl) + M{\p\) ' 

The renormalization function Z^(|p|), the mass function 
M(|p|) and the "energy form factor" a(|p|) are derived 
in Appendix C: 

-I -1 



ki\P\. 



a(lpl) = 



M{\p\) - 



dpi 
27r 



M\p\,Pi) 



p dp4_ D 



2ir 



]P|,P4) 



dp4 4 



]P|>P4) 



(25) 



where the integrals here simply indicate a statistical av- 
erage over the energy, given the p4 independence of the 
dressing functions ^'s and B discussed in Sec. IV B. 

To check renormalizability we now need to establish 
the scale invariance of M(|p|), a(|p|) and study the 
scaling properties of Z(;-(|p|). Notice that M and 
can also be extracted from the equal-times propaga- 
tor (cf. Eq. (14)), while a always requires the full p^- 
dependent propagator (cf. Eq. (24)). 



In Fig. 3a we show the mass function Af(|p|) from 
configuration sets (a)-(d), with scale a sa 0.12 fm, com- 
pared to configuration sets (/) and (g), which have a 
scale a = 0.086 fm. These sets are chosen to have ap- 
proximately the same physical volume, so as to minimize 
finite size effects at the two different cutoffs. As can be 
seen M(|p|) nicely agrees for the sets with similar masses, 
namely (a) (15.7 MeV) and (/) (14.0 MeV) on one side 
and (6) (31.5 MeV) and (g) (27.1 MeV) on the other side. 
We compare thus in Fig. 4b the corresponding wavefunc- 
tion renormalization functions Z^(|p|) from configuration 
sets (g) and (6), finding a good agreement once we nor- 
malize both to -^^(C) = 1 for the scale C, = 3.0 GeV. 
The scaling of the Landau gauge quark propagator was 
checked on the same lattices in Ref. [32], giving very 
similar results. 

In Figs. 3b and 4c we compare M(|p|) and Z^(|p|) from 
configuration sets (/) (14.0 MeV) and (i) (11.9 MeV). Al- 
though on the latter we have less statistics and could only 
calculate Kogut-Susskind fermions, the agreement is still 
quite good. 

We can thus conclude that the static propagator 
Eq. (14) is multiplicatively renormalizable. It should be 
stressed, however, that the wave function normalization 
(and a, cf. next paragraph) show satisfactory scaling 
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(b) M from (/) and (i) 
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fm, m = 14.0 MeV. dynamical 
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from (/) and (h) 

FIG. 3. Scaling of the mass function M 
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|k| [GeV) 

(d) M from 



only when the improved Asqtad action is taken. This is 
in contrast to the running mass M(|p|), which is very 
robust against lattice artifacts, cf. e.g. Fig. 3d, where the 
unimproved Wilson action already gives a reasonable re- 
sult (m = 212 MeV for all three data sets). A similar 
effect was also noticed in Landau gauge [32, 34-37]. 

Turning now to the full propagator Eq. (18), Figs. 5b- 
5c show the function a, i.e. the ratio of At to As, for 
different configuration sets. The scaling behaviour is very 
good both for the dynamical configurations Fig. 5b and 
for the quenched configurations Fig. 5c. We can thus 
conclude that the full (p4-dependent) propagator Eq. (18) 
is also multiplicatively renormalizable. As argued above, 
this fact allows us to define a dispersion relation for a 
single pseudo-quark as in Eq. (19). 

Interestingly, Zi^ turns out to be much more suppressed 
in the IR than its Landau gauge counterpart, cf. Fig. 6b, 
while a has a mild momentum dependence. The result- 
ing dispersion relation Wf^dpj) is then found to be IR 
enhanced, as can be seen in Fig. 5d for sets (d)-(e), for 
which we have the best signal to noise ratio. If Z'^/a 
should indeed prove to vanish in the IR according to a 
power law, 

Zl{\p\)/a{\p\) (X 

we would find for the quark dispersion relation a be- 
haviour similar to the Gribov formula for the gluon 



[7, 13, 14]: in the ultraviolet ujf oc \p\ and in the infra- 
read a power law utp oc IpI"*^. The last relation describes 
an infrared diverging effective (pseudo-) energy for the 
quark, which would extend the Gribov-Zwanziger con- 
finement scenario to the fermionic sector of QCD. With 
our present data we find a behaviour compatible with an 
infrared exponent of fc « 0.25, but of course a reliable 
estimate would need better statistics on larger volumes, 
which we plan to collect in the near future. Notice also 
that the minimum of Wi?(|p|) is around |fc| ~ 1.2 GeV, 
which is compatible with the at ^ extrapolation of the 
gluonic Gribov mass [16] and the estimate of A we will 
extract in Sec. IV F. 



D. Effects of dynamical quarks 

The dynamical and quenched configuration sets (/) 
and (h) have same size and scale. This makes them well 
suited to study quantitatively the effects of quenching. 
To this end we have calculated on both sets the Coulomb 
gauge quark propagator with equal bare mass, viz. the 
light sea quark mass of set (/). The resulting mass and 
renormalization functions M(|p|) and Z^(|p|) are com- 
pared in Figs. 3c and 4d. We find a slight screening of 
the dynamical mass generation in the case of full QCD, 




m = 14.0 MeV, IPG > 
m . 27.1 MeV, IPG > 
m . 15.7 MeV, IPG > 
m . 31.5 MeV, IPG > 
m = 47.3 MeV, IPG > 
m » 63.1 MeV, IPG ' 

2.5 3 3 



|k| [GeV] 




(a) Z from (a)-(g). 



(b) Z from (b) and (g) 




64'x144, a=0.06 fm, m = 1 1 .8 MeV • 
28^x96, a = 0.086 fm, m = 14.0 MeV > 

3 4 5 

|k| [GeV] 



(c) Z from {/) and (i) 




28 x96, a-0.086 fm, m - 14.0 MeV, dynamical ^ 
28^x96, a-0.086 fm, m = 14.0 MeV, quencfied ' 

2 3 4 

|k| [GeV] 

(d) Z from {f)-(h). 



FIG. 4. Scaling function Z. The renormalization point is set at 3 GeV in (4a) and (4b), 3.7 GeV in 4c and 4.64 GeV in 4d 



similarly to what was seen in Landau gauge [37] 



E. Comparison to Landau gauge 

In Fig. 6a and 6b we show the comparison between the 
mass and renormalization function in Coulomb gauge, 
Eq. (14), and Landau gauge, Eq. (11). The data for the 
Landau case was taken from Ref. [32]. Up to slight dif- 
ferences in the intermediate and high momentum region, 
which are however still within our error bars, M is almost 
identical in both gauges, especially in the IR, while 
shows a much stronger gauge dependence. For this com- 
parison, the renormalization point has been set at C = 
4.64 GeV. These findings nicely agree with the discussion 
in Sec. IIIB. 



F. Chiral limit 

For the configuration sets (a)-(d) the dynamical 
masses decrease at constant cutoff. We have therefore 
attempted a chiral extrapolation of the mass function 
determined for these sets with Asqtad fermions at the dy- 
namical point. This procedure should minimize system- 
atic errors due to partial quenching. The result is shown 



in Fig. 6c. In addition, we have also made a chiral ex- 
trapolation of the set (i) using Kogut-Susskind fermions 
(Fig. 6d), where the mass has been fixed to five different 
values, with only the last at the dynamical point. The re- 
sult of the latter, although quite noisy in the IR, allows 
us to extend the momentum range and extract further 
information in the UV. Within error bars the two chiral 
extrapolations agree, cf. Fig. 7a. We find that all data 
can be well described by a function of the form: 



M{\k\,mb) = 



,{mb) 



log (^e 
mr{mii) 



A2 



log (^e 



A2 



7 



(26) 



where 6 is a constant, 7 is the anomalous dimension and 
A is the QCD scale. In this fit we have also defined 
a renornialized Coulomb dressed mass and a chiral 
quark mass, both of which depend on the bare mass mf,. 
Performing a simultaneous fit for all data we find b ~ 
2.9(1), 7 = 0.84(2) and A = 1.22(6) GeV with xVd-o.f.= 
1.06. For the extrapolated data of Fig. 7a we further 
find (constraining = 0) = 0.31(1) GeV, in very 
good agreement with the constituent mass expected from 
quark models. In Fig. 7b we show the dependence of 
on the bare mass mi,. For bare masses larger than mb — 



9 



1 

0.5 


-0.5 
-1 



A,(p ) Coulomb gauge 
A((p^) Coulomb gauge + IPG 



0.5 1 1.5 2 2.5 3 3.5 4 
|k| [GeV] 

(a) At from (a) 

28^x96, a = 0.086 fm, m = 14.0 MeV, IPG, quenched □ 
20^x64, a = 0.121 1m, m = 63.0 MeV, IPG, quenched — »^ 



|k| [GeV] 
(c) a from (e), {h) 



0.8 
0.7 
0.6 
0.6 
= 0.4 
0.3 



28''x96. a . 0.086 Im, m . 27.1 MeV, IPG > 

20^x64, a = 0.121 1m, m = 31.5 MeV, IPG > 

20'x64, a = 0.120 1m, m = 47.3 MeV, IPG > 

20'x64, a = 0.1 19 fm, m = 63.1 MeV, IPG > 




I 1.5 2 2.5 3 

|k| [GeV] 

(b) a from {a)-{h) 



quenched 
o unquenched 




0.5 



|k| (GeV) 

(d) ujp from {d)-(e) 



FIG. 5. Scaling of the temporal function a; ujf as in Eq. (19). 



100 MeV the chiral quark mass practically disappears; on 
the other hand we find that the Coulomb gauge dressed 
mass nir increases more than linearly with mf,, cf. Fig 7c. 
One should not however over- interpret such result, since 
the splitting between the two masses in the intermediate 
region will strongly depend on the fit function chosen. 
A more reliable extraction of will need to be based 
on data at smaller lattice spacing, pushing the fits in 
the higher UV region. Finally, Fig. 7d displays the total 
constituent IR mass M(0) = + from Eq.(26) as a 
function of rrifc; this should turn out to be fit independent 
as soon as the available data will push far enough in the 
IR region. 

Similar calculations with Asqtad fermions at the finest 
spaced lattice and with the largest overall volume are 
currently under way to confirm our results. 



counterpart. However, for future work, it will allow for 
a more efficient determination of the chiral parameters, 
the constituent quark masses and the quark anomalous 
dimension, since as we have shown the Dirac operator in 
Coulomb gauge is energy independent and its static ver- 
sion, which holds all information about the running mass 
and the wave function renormalization, can be inverted 
at each fixed time slice. 

Moreover, the full (non-static) Coulomb propagator 
gives access to an effective energy for the quark field. 
With our current data, the dispersion relation is compat- 
ible with a Gribov kind of behaviour. Quark confinement 
could thus be explained in this picture by a diverging IR 
effective energy, which is a first evidence that the Gribov- 
Zwanziger confinement scenario could be extended to the 
fermionic sector of QCD. 
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(d) R4 from (i) 



FIG. 6. Landau vs. Coulomb gauge (top); chiral behaviour of Af (bottom) 



Appendix A: Fixing the residual gauge 

After fixing the Coulomb gauge at each time-slice X4 we 
are still left with the freedom to perform a global gauge 
transformation 3(0:4) at each X4. To determine it we im- 
plement the lattice equivalent to the residual continuum 
condition (7). Let us begin by averaging, for every fixed 
time-slice x^, over the spatial coordinates x in 1/4 (x), 



u{x4) 



(Al) 



Next we project the matrix u(x4) back to SU(3) by 
Cabibbo-Marinari cooling [53], i.e. we look for a u{x4) G 
SU{3) satisfying: 



maxCHe tr [m(x4)u^(x4)] 



(A2) 



Although such projection is not unique, we believe it 
to be a "natural" choice; of course, different projections 
will lead to different temporal gauges. However, since all 
dressing functions in S turn out to be energy indepen- 
dent, the choice of a different prescription should even- 
tually not influence the results. 

We now seek a gauge transformation g(x4) such that 



u{x4) — > const. 



(A3) 



This can be achieved defining the integrated Polyakov 
loop: 



P = tr 



J|u(a;4) 



and choosing g{x4) such that for all 2:4 

g{x4)u{x4)gHx4 + l) = P'^^, 



(A4) 



(A5) 



where T is the temporal extension of the lattice and P^^^ 
is the T-th root of P. To do so we (arbitrarily) choose 
g{0) = 1 and then determine all g{x4) recursively, 



^\x4 + l)^uHx4)gHx4)P'^^- 



(A6) 



If we now gauge all links via the transformation g{x4) we 
have 



U^{X,X4) 
U-{X,X4) 



g{x4)U4{x,X4)g^X4 + 1) (A7) 
g{x4)U,{x,X4)g^{x4), (A8) 



respectively. The new temporal links U!j^ obey Eq. (A3); 
if u ~ u as in the case of SU (2) , this condition would via 
Eq. (Al) translate into 84 J2x ^^4(^)1 which is the equiv- 
alent to the continuum condition (7). For G = SU{3), 
however, the sum of colour matrices is not proportional to 
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FIG. 7. The chiral limit taken from the dynamical point of configuration sets (o)-(d), compared to the chiral limit from set (i) 
(partially quenched) and to the fit to Eq. (26). 



an SU(3) matrix, and we can only choose u(x4) as close 
as possible to u{x4). The cooling step Eq. (A2) imple- 
ments thus the continuum condition Eq. (7) as much as 
possible for the color group G = SU{3). From Eq. (A8), 
we also see that the new transformation g{x4) acts as a 
global gauge transformation at each fixed time slice X4, 
i.e. the Coulomb gauge condition remains unaffected. 



Appendix B: Staggered quark propagators 

We extend here to Coulomb gauge the decomposition 
of the staggered quark propagator in Landau gauge per- 
formed in [37, 54] . The Kogut-Susskind propagator reads 
at tree-level: 

^a^iiy^ = + moSaf} (Bl) 



where a, /3 are staggered multi-indices, a — (ai, . . . , 04), 



e {0, 1}, while = \[ 5^ 



and 



(7;.)a/3 = (-l)"^'5a+0(M),/3 



(B2) 



where the v*"^ component of the 4- vectors 9{fi) is given 
by: 



1 if v < fi, 



else , 

giving the staggered Dirac algebra 



(B3) 



(B4) 



The discrete momenta 5^ = 0(7^ are restricted for stag- 
gered fermions to the inner half of the Brillouin zone, 
Qfj, £ (~fjf]i ^^'^ related to the ordinary lattice 
momenta p^t G (— 7r,7r] via 



(B5) 



In the following we will use the common abbreviation 
= sinq^. 

Summing over one of the propagator's multi-indices, 
E/3 = Ek....,/34=0' define 



(B6) 
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In order to evaluate the r.h.s. of Eq. (B6), we note that 

We now want to extract the dressing functions of the 
inverse Coulomb propagator: 

G-\q) = i{-irhaAs{\q\,ql)+ti-ir-k4aAt{\q\,ql) 
+ i(-l)".+"4iaAd(|q|, «|) + BM,<ll) • (B8) 

On the other hand, we can also express the propagator 
in terms of new dressing functions As, At, Ad and Bm- 

Ga{q) ^ ~t{-irkaAs{\q\,ql) - i{-ir'kaAt{\q\,ql) 

- i{-ir+'''kaAdi\q\,ql) + BM,(ll) 

= [ _ t{-l)"^kaAA\q\,ql) - ti-l)"' ha At{\q\, qj) 

~z{-ir+'^-kaAdi\q\,ql)+BMqi)]/D^{\q\,ql) 

(B9) 

where we have defined 

D'{\q\,qi)^Y.''i^''^'{\Q\,Ql) + kWAU\q\,ql) 

i 

+ J2kyAl{\q\,ql) + Bli\q\,ql). (BIO) 



After some algebra it can be shown that 
D\\q\,ql)V\\q\,ql) = l 

where 



(Bll) 



V\\q\,ql) ^ J2 kfa^AU\q\,ql) + kja^AUlqlql) 

+ Y,kfa'Al{\q\,ql)+Bl{\q\,ql). (B12) 

i 

Now we multiply Eq. (B9) by (-1)"', (-1)"", 
or 1, respectively, and sum over a using 



(B13) 



Taking finally the trace with respect to color indices we 
obtain: 



DH\q\,ql) 16iV,E,fc?« 

xEE(-ir4itr[G„(g)], (B14) 



At{\q\,qi) 



2x _ At{\q\,ql) _ i 



xE(-l)"44tr[G„(5)], (B15) 



X J2 , (B16) 



(B17) 

From Eq. (Bll) and Eqs. (B14-B17) we can now easily 
extract the dressing functions Ag, At, A^ and B^- 

Turning next to Asqtad improved fermions, their tree- 
level propagator reads: 



1 2/ ^ 

1 + -sm [qf,] 



■ rh5„p (B18) 



and the decomposition given above can be performed in 
essentially the same way. We only have to keep in mind, 
though, that the dressing ftmctions will also get contri- 
butions from the tadpole factors uq, which have to be 
eliminated a posteriori. 



Appendix C: Formulas for Z, M and a 

To extract 2'(|p|), M(|p|) and a we can proceed as 
in Appendix B. Recalling that ^4^ = we write the 

propagator as (fc^ = sinp^): 

S'-i(p) = i^As{\p\,pi) + i^iaAt{\p\,p^) + S„(|p|,p4) 

(CI) 

As{\p\,P4) , .„ At{\p\,PA) , B„,i\p\,P4) 



'^v^qpm + v^\p\,p4) 



(C2) 



and 



Sip) = -i^As{\p\,P4) - i^4aAt{\p\,P4,) + Bm.i\p\,P4) 

(C3) 

- _ -M^ jMWlM _ ill. „ M\P\,P4) , Bm{\p\,P4) 

- '^D^{\p\,P4) '''"'d^{\p\,P4)^ D^\p\,p4)- 

(C4) 

Here the denominators are 

D^\p\,P4) ^ k'a'Al(\p\,p4) + kyAl{\p\,p4) 

+ BI{\pIp4), (C5) 
VWp\,p4) = k'a^Ali\p\,P4) + kyA^ti\PlPi) 

+ Bl{\p\,p4), (C6) 

and they still obey 

D\\p\,pi)V\\p\,pi) = l. (C7) 
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This looks just like the Landau gauge case, but with 
an extra structure function. Applying the same line of 
reasoning as in Ref. [37, 54] now leads immediately to 
Eq. (25) in the main text. 

Turning now to the static propagators, we first note 
that At(|p|,p4) is an even function of P4, so that the 
74 contribution cancels by T-symmetry when integrating 
S~^{p) over the (energy) Brillouin zone. Therefore: 



s-\p)=i^ r '^M\p\,P4)+ f ^B™(ipi,P4) 



(C8) 



J —IT 



dp4 As{\p\,P4) , r dpi Bm{\p\,P4) 



and 



l3m{\p\,P4) 



= — i^a J 



(CIO) 



dp4 As{\p\,p4) r dp4 B^{\p\,pi) 



27r D^{\p\,p4) 



+ 



f 

J — TT 



27r D^{\pIp4) ■ 

(Cll) 



From Eq. (C8) we obtain: 



S{P) 



-i1feaS\ ^AMp4) + ^B™(|p|,p4) 



p'\\p\) ^ ea^ 



' dp4^ 
n 27r 

277 



p'{\p\: 



As{\p\,P4) 



Bm{\p\,P4) 



(C12) 



(C13) 



while inverting Eq. (CIO) yields: 

f^A(|p|,P4) + 1% ^BM,P4) 



s-\p) 



V^\p\) 



V\\p\)^k\' (^J^^^As{\p\,P4) 



" dp4 

277 



^^(1^1,^4) 



Prom Eq. (C8) and Eq. (C14) we thus obtain: 



^ dP4 „ , . /:.^gm(iPi,P4) 

^BM,P4) = 



while Eq. (CIO) and Eq. (C12) give: 

dP4 , , I^'^M\PIP4) 
^Bmi\plP4) = . 

Writing the static propagator as: 

-ilta + Mi\p\) 



S{p) = Z{\p\[ 



fe2a2+M2(|p|)' 
S-'{p) = Z-\\p\){iiki + M{\p\)). 



(C14) 
(C15) 

(C16) 
(C17) 

(CIS) 
(C19) 

(C20) 
(C21) 



and comparing Eq. (C14) and Eq. (C21) we finally have: 



M{\p\) 



/-.frA(|p|,P4) i:^tA{\p\,P4) 



(C22) 



i:.tt^m{\PiP4) i:^^B^{\p\,p4 



i:^^M\p\,P4) j:^tA{\p\,P4) ■ 

(C23) 
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